Numerical signs for a transition in the 2d Random Field Ising 

Model at T = 



Carlos Frontera and Eduard Vives 
Departament d'Estructura i Constituents de la Materia, Facultat de Fisica, 
Universitat de Barcelona, Diagonal 647, E-08028 Barcelona, Catalonia, Spain. 

and 

Escola Universitaria Politecnica de Matard 
Av. Puig i Cadafalch 101-111, E-08303 Matard, Catalonia, Spain. 

(February 1, 2008) 



Abstract 



Intensive numerical studies of exact ground states of the 2-d ferromagnetic 
random field Ising model at T = with gaussian distribution of fields are pre- 
sented. Standard finite size scaling analysis of the data suggests the existence 
of a transition at a c = 0.64±0.08. Results are compared with existing theories 
and with the study of metastable avalanches in the same model. 
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The study of systems with quenched disorder has been a challenging problem since the 
last 15 years. The interplay between thermal fluctuations and disorder has a great influence 
on the existing phase transitions. Many systems are known to exhibit such phase diagrams 
highly determined by the degree of disorder (vacancies, impurities, dislocations, etc.) Among 
other, the most typical examples can be found in magnetism, superconductivity, structural 
phase transitions, etc. For such systems different models have been proposed. The Ising 
model with quenched disorder is one of the simplests and it has the advantage that the pure 
model is well known. The disorder can be of two types: (i) symmetry breaking terms like 
random-fields or random magnetic impurities, and (ii) non-symmetry breaking like random- 
bonds, vacancies, etc. For all the cases different probability distributions of disorder have 
been studied. Here we will focus on the study of the Random Field Ising Model (RFIM) 
in two dimensions (2d) with a gaussian distribution of fields. Since many years ago there 
has been a discussion concerning the possibility that, for such a 2d model with symmetry 
breaking random fields, it exists or not order at low temperatures. The initial studies lead to 
a certain controversy: the Imry-Ma |IJ argument suggestes that the lower critical dimension, 
below which ferromagnetic order is distroyed, is di < 2 with d = 2 being the limiting case. 
Renormalization group expansions @ around d = 6 lead to the "dimensional reduction" 
argument suggesting that di — 3, discarding the possibility for ordering in the 2d- RFIM. It 
has also been suggested || that there are new types of order for d > 1. This controversy is 
probably due to the difficulty in balancing the two ingredients of such models: disorder and 
thermal fluctuations. 

More recently a new approach to disordered systems has been proposed, namely the 
study of disordered systems at T = 0, i.e. withouth thermal fluctuations. From a theoretical 
point of view this simplifies the problem withouth making it trivial. Moreover, several 
experimental systems exhibit phase transitions that can be catalogued into this "athermal" 
category: two examples are ferromagnetism at low temperature under an external magnetic 
field M, and martensitic transformations ||. Both systems present a first-order phase 
transition that can be crossed by sweeping a control parameter and are greatly affected 
by the presence of quenched disorder. We will concentrate on the study of the 2-d RFIM at 
T = for different values of the standard deviation a of the gaussian distribution of fields. 
Our goal is to look for signs for the existence of a phase transition at a certain a c from a 
ferromagnetic ordered state for o < o c to a disordered state for a > a c . For the 3d-RFIM 
at T = 0, ground state studies ||[7j and renormalization group arguments || reveal the 
existence of such phase transition, but to our knowledge no results for the 2d case have 
been published. Figure [1] summarizes the finite size scaling study presented in this paper. 
Data corresponds to estimations of a C L obtained using different methods, as a function of 
the linear system size L 1 ^ (where 1/u — 0.5 is the exponent characterizing the correlation 
length divergence). The standard extrapolation to L — > oo, as will be discussed, renders 
a c = 0.64 ± 0.08 different from zero. 

We consider the 2d-RFIM on a L x L square lattice with periodic boundary conditions 
and with the hamiltonian H = — Yn'J' ^i^j ~ ££i Sih h where i and j are indices sweeping 
the full lattice = 1, . . . , N = L x L), the sum refers to nearest-neighbours (n.n.) pairs, 
Si = ±1 are spin variables and hi are independent random fields distributed according to 
the gaussian probability density with (h) = and (h 2 ) = a 2 . The advantage of using a 
continuous distribution is that, for almost any configuration of fields [hi] the ground-sate 
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is not degenerated. The order parameter is the magnetization of the system defined as 
m = J2Si/N. Because the ground state is unique, thermal averages are meaningless. Since 
we are interested in the dependence of the system properties with the amount of disorder 
a the only possible averages with physical meaning are the ensemble averages {(■■■)(&)) 
performed over different realizations of the random fields with a certain fixed degree of 
disorder a. Experimentally this has to be understood as averaging measurements on different 
samples which have been prepared with the same amount of disorder. 

The zeroth-order mean field (MF) theory was solved 15 years ago |J. A solution with 
the order parameter (m)(cr) 7^ appears for a < o c = -^=, and the phase transition is 
continuous. Of course this MF result cannot be expected to be correct, neither to reflect 
any dependence on dimensionality. Moreover, the MF studies can be extended to higher 
orders by exactly treating larger and larger clusters of spins. For thermal phase transitions 
this is known to extrapolate to the exact value of the critical temperature. The first-order 
approximation is the Bethe approximation which consists in solving exactly a cluster of a 
central spin and its four n.n. The method can be extended to larger clusters. We have found 
a continuous phase transition at: cr c = 2.76, 2.48, 2.12 and 1.98 for clusters of N = 5, 13, 25 
and 41 spins respectively. These results are indicated with stars in Fig. [I] (considering 
L = y/N). 

A better approach consists in looking for exact ground states by using the max-flow 
min-cut theorem || . We have designed an algorithm which solves a set of different a values 
with a minimization time that grows like L 4 . We have studied lattices with L = 4,8,16,32,64 
and 128 and taken averages over 10 5 , 10 4 , 10 4 ,5 10 3 , 10 3 and 30 realizations of random fields 
respectivelly. The inset in Fig. [I] shows a typical example of a ground state for a = 1.0 and 
L = 64. We have focused on the computation of different magnitudes. The order parameter 
has been estimated from (|w|)i(cr) and J (m 2 ) L (a). The subscript L indicates that such 
quantities will, in general, depend on system size. We have also measured the susceptibility 
as Xl(ct) = N {(rn 2 ) L — (|m|)^J (which, for large a, tends to 1 independently of L) and the 

Binder's cumulant <7l(c) — 1 — ( m4 ) lI (3( m2 ) 1) ■ Also the correlation length £l(ct) can be 
computed by fitting an exponential decay to the spin-spin correlation function. 



Figure |] shows the behaviour of y (m 2 )^ as a function of a for different system sizes. One 
can estimate <j c l as the inflection point of a fitted third order polinomium. Data can be scaled 

where M 



using the standard finite size scaling assumption: y ( m2 ) l ~ L@M L l l v (o — <j c l) 
is the corresponding scaling function. The exponents j3 and v can be estimated by fitting the 



power-laws yj \m 2 )(a = a cL ) ~ L 13 and — s — (a = a cL ) ~ L p+1/v . One gets (3 = -0.038 ± 

0.009 and \ jv — 0.54±0.04. The scaled data is shown in the inset in Fig. |2|. A similar scheme 
can be applied to the study of (|m|)i,(cr) rendering j3 = — 0.026±0.017 and \ jv — 0.54±0.05. 
Susceptibility Xhip)-, shown in Fig. ^, exhibits a peak at o = <j c l which shifts and increases 
when increasing L. Data can also be scaled using \l ~ L a x L x l v {o — a c i) . Power law fits 
to the heigth and curvature of the peak render a = 1.89 ±0.03 and \ jv = 0.46 ±0.05. Scaled 
data is shown in the inset of Fig. |3|. Figure ^ shows the behaviour of £l(ct). The peak gives 
an independent measure of 0&,. The continous line is an estimation of £(cr) for L — > 00 that 
will be discused later. Note that the behaviour of the curves is compatible with the finite 
size scaling hypothesis, i.e., £^ follows the behaviour corresponding to the infinite system 
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up to a certain ^ max = KL. (Data is compatible with K ~ 0.08). A last estimation of a c i 
can be obtained from ^l(ct). For all the studied sizes <?l(ct) takes a low value for large a and 
reaches the value 2/3 at a certain a C L. This estimation is independent of L as suggested in 
ref. 0. We want to note that (3 ~ 0, which means that the order parameter increases very 
fast to 77i — 1 after the transition. A low-cr first-order expansion renders 1 — m ~ 10~ n for 
a = 0.6. This fact may also explain why it is very difficult to measure with numerical 

accuracy enough to check for a crossing point which is the standard procedure to locate the 
transition. Note also that (3 ~ and a ~ 2 would suggests that m exhibits a lack of self 
averaging [jTIJ. 

The different estimations, as explained in the previous paragraph, of are plotted in 
Fig. in front of L~ x l v . The open symbols correspond to the estimation from (\m\)L{o') 



(m 2 )^ (□), Xl(&) (O), (A) and <?l(o") (y). 111 order to extrapolate the data to 

L — > oo, we have used the standard expansion for the divergence of £, up to second order: 
£ ~ (a — <7 C )~ U [1 + C(cr — <j c )}. Now, imposing that a C L is determined by the condition 
£ = f^L one gets: a cL = a c + C X L~ X I V + C 2 L-^ U . Such parabolic fits are also shown in Fie;. 
[I] with continuous lines. The extrapolated a c lay all within a c = 0.64 ± 0.08. To get an idea 
of the error margins, we have also fitted the first order expansion (C 2 = 0) leaving v free, 
rendering a c = 0.65 ±0.1 and v = 1.8 ± 0.2; or fixing v = 2 which renders cr c = 0.56 ± 0.06. 

The existence of this phase transition is in apparent contradiction with previous results. 
It has been proved |Tl|] that the RFIM has a unique Gibbs state in the thermodynamic limit, 
i.e. for a given configuration of the random fields the ground state is unique. This can be 
misunderstood [0 as a proof that the ordered phase cannot exist. When considering the 
ensemble of all possible realizations of the random fields corresponding to a certain value 
of a it may well be that the distribution of magnetizations changes from a single peak one 
(for large a) to a bimodal one for small values of o. Thus the phase transition we are 
proposing should be understood as existing in this ensemble rather than for a single system 
for which the ground state is unique. It is true that there is an open question here concerning 
the size of this ensemble: in the thermodynamic limit, is there more than one realization 
of disorder compatible with a certain a? We understand that given the discreteness of the 
Ising lattice and the continuity of the random fields one can still consider the existence of 
such an ensemble. 

Another interesting point is the comparison of our results with studies suggesting an 
exponential divergence of the correlation length £ ~ exp{— B/(a — <J C ) T }. On the basis of 
the study of the interfaces separating regions with m > and m < Binder |13[ derived a 
theory with r = 2 and a c = 0. We have tested the validity of this theory by studying the 
corresponding finite size scaling hypothesis [r5,l4]|: l/cr? L oc ln(L~ l ). The inset in Fig. |] 
shows the comparison of this behaviour with the standard one we propose in this work on 



top of the data corresponding to the estimations of <t c l from the inflection point in y (m 2 )i 
(both fits have two free parameters). The standard theory works better. We have also 

L l / y + A)~"), the finite size 



tested that, assuming the standard hypothesis (£l ~ 
scaling of £l is better than using Binder's hypothesis. Moreover, Binder's theory proposes 
that for large enough systems the configurations with total m ~ will be more and more 
frequent. We have not observed the existence of many of such "slab" configurations but 
found ground-states with closed domains like that in the inset of Fig. |]. Figure |5] shows the 
probability P(m) obtained from the computations of a very large number of exact ground- 
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states for a system with L = 64. Clearly, for a < er c , the configurations with m ~ have 
much less probability than the configurations with m ~ 1. The reason for the failure of 
Binder's theory could be that, in order to perform the thermodynamic limit he uses a very 
anisotropic system with open boundary conditions. 

Our data can be compared with the studies of the evolution of the RFIM at T = 0, 
obbeying a local relaxation dynamics. It has been found that when sweeping the external 
field the system evolves by avalanches between metastable states. At a certain degree of 
disorder o c the distribution of avalanches becomes critical. In Fig. [I] we show the values of 
the <t c l corresponding to the 2d case from Ref. The behaviour is very similar to the 

equilibirum data. Different extrapolations to L — > oo have been reported (cr c = 0.75 ± 0.03 
|i~5fl , er c = 0.54±0.04 fl6|]) but all are close to the equilibrium one. Concerning the exponents, 



for the metastable studies the exponent (3 has also been found to be very small, while previous 
reported values for v are 1.6 ± 0.1 [15] and 5.3 ± 1.4 |16|]. Therefore we suggest that the 



metastability phenomena found in the out-of-equilibrium studies might be associated with 
a real underlying equilibrium phase transition at o < o c for zero external field. It should 
be mentioned that in the context of these out of equilibrium phase transitions Sethna and 
collaborators have |i| proposed a theory with exponential divergence of £ with r = 1 and 
a c = 0.42 ± 0.04. Our data is not consistent with such theory. If we perform a fit in the 
evolution of <j c l{L) leaving a c and r free we get r = 0.6 and a c = 0.25. We can obtain still 
a good fit (and good scalings) by taking r = 1 and a c = 0, although we cannot provide 
any physical explanation for such a behaviour. We finally want to point out that the phase 
transition we have found at T = may also be related to the change in the type of growth 
found at o = 0.33 in the studies of the depinning transition in the same model [[17] . 



In conclusion, we have presented a finite size scaling analysis of numerical data for systems 
up to L = 128, that suggests that the RFIM with a gaussian distribution of fields, at T = 
exhibits, a phase transition at a c = 0.64 ± 0.08. The ensemble average of the magnetization 
changes from (m) = for a > a c to a state with (m) ^ for a < a c . The transition is 
characterized by the exponents l/v = 0.5 ± 0.05, (3 = -0.03 ± 0.02 and a = 1.89 ± 0.03. 
The possibility of a exponential divergence £ ~ exp(5/a) cannot be excluded. 

We acknowledge finantial support from CICyT project number (mat95-0504), CRAY 
and FCR. OF. also acknowledges the Comissionat per a Universitats i Recerca (Generalitat 
de Catalunya) for finantial support. 



5 



REFERENCES 



Y. Imry and S.K. Ma, Phys. Rev. Lett. 35, 1399 (1975). 

See G.Parisi and N.Sourlas, Phys. Rev. Lett. 43, 744 (1979) and J.Bricmont and 

A.Kupiainen Phys. Rev. Lett. 59 (1987). 

E.B.Kolomeiskii, JETP Letters 52, 538 (1990). 

P.J.Cote and L.V.Meisel, Phys. Rev. Lett. 67, 1334 (1991). 

E.Vives, J.Ortm, L. Manosa, I. Rafols, R. Perez-Magrane and A. Planes, Phys. Rev. Lett. 
72, 1694 (1994). 

A.T. Ogielski, Phys. Rev. Lett. 57, 1251 (1986). 

M.R. Swift, A. J. Bray, A.Maritan, M.Cieplak and J.R.Banavar, Europhys. Lett. 38, 273 
(1997). 

M.E.Newman, B.W.Roberts, G.T.Barkema and J.P.Sethna, Phys. Rev. B 48, 16533 
(1993). 

A. Aharony, Phys. Rev. B 18, 3318 (1978). 

A.M.Ferrenberg, D.P.Landau and K. Binder, Journ. of Stat. Phys. 63, 867 (1991). 
M.Aizenman and J.Wehr, Phys. Rev. Lett. 62, 2503 (1989). 

T.Nattermann, "Theory of the Random Field Ising Model", to appear in Spin Glasses 
and Random Fields, ed. P.Young, World Scientific (see |cond-mat / 97052951) . 



K.Binder, Z.Phys. B 50, 343 (1983). He studies a RFIM with the ±h distribution but we 
think that his arguments can also be applied to the present gaussian distribution. 
The same scaling has been proposed in: A. J. Bray and M.A.Moore, J. Phys. C 18, L927 
(1985). 

E. Vives, J. Goicoechea, J. Ortm and A. Planes, Phys. Rev. E 52, R5 (1995). 
O.Perkovic, K.A.Dahmen and J.P.Sethna, to be published (1997). 
Hong Ji and M.O. Robbins, Phys. Rev. A 44, 2538 (1991). 



6 



FIGURES 




FIG. 1. o c l versus L~ x l v . The results have been obtained using MF at zero and higher orders 
(stars), exact solution of finite lattices (o, □, O, A, v) an d studies of the metastability behaviour 
(black triangles) from Ref.[15]. Typical error bars are displayed. The inset shows an example of 
ground-state of a L = 64 system with a = 1.0. 
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FIG. 3. Behaviour of Xl(<t) for L = 16,32,64 and 128. The inset shows the same data scaled 
using a = 1.89 and \ jv = 0.46. 
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FIG. 4. Correlation length £l(<t) for systems with L = 16,32,64,128. The continuous line 
shows the behaviour of £ = ^4(u — <T c ) _Jy [l + C(cr — a c )\ with a c = 0.64 and z> = 2. The inset 
shows the finite size dependence of a~1 versus /n(L _1 ). Data is the same as in Fig. 1 (□). The 
continuous line is the standard scaling used in this work and the dashed one is the best fit of the 
theory by Binder et al. in Ref.[13]. 
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FIG. 5. Histograms of the magnetization distribution for a L = 64 system and for 
a = 1.20,1.15,1.10,1.05 and 1.00 (from bottom to top). The distribution evolves from a sin- 
gle peaked one at a > a c to a two-peak distribution for a < a c . The curves are shifted 0.2 units 
each to clarify the plot. 
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